package scu.maqiang.cfd;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2Q41;
import scu.maqiang.mesh.Mesh2Q4;
import scu.maqiang.mesh.MeshFormat;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;
import scu.maqiang.possion.Possion2DQ41;

public class PotentialFlow {

	public static void main(String[] args) {
		String meshFile = "mesh2D\\quarterCircle.msh";
		Mesh2Q4 mesh = new Mesh2Q4();
		mesh.readMesh(MeshFormat.msh, meshFile);
		mesh.toTecplot("FlowProblem\\quarterCircle.dat");
		
		Possion2DQ41 psQ4 = new Possion2DQ41(mesh);
		psQ4.GeneralPossionProblem(BVPType.COMMON, 1.0, 0.0, new int[] {3}, 40.0,
				new int[] {5}, -5.0, null, 0.0, 0.0, "FlowProblem\\potentialFlow.dat");
		
		FES2Q41 fs = new FES2Q41(mesh);
		SRMatrix A = new SRMatrix(fs.GetNdof());
		double[] RHS = new double[fs.GetNdof()];
		fs.assembleStiff(new double[] {1.0}, BVPType.COMMON, A);
		fs.applyBC_MBN(A, Direct.X, 2, 4);
		fs.applyBC_MBN(RHS, Direct.X, 0.0, 2);
		fs.applyBC_MBN(RHS, Direct.X, 40.0, 4);
		NewIterSSolver solver = new NewIterSSolver(A);
		double[] x0 = new double[fs.GetNdof()];
		solver.PCGSSOR(RHS, x0, 1.5, 1);
		double[][] gradU = new double[2][fs.GetNdof()];
		fs.computeGradient(x0, gradU);
		System.out.println(MVO.toString(gradU));
		mesh.toTecplot("FlowProblem\\StreamFlow.dat", gradU);
	}

}
